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ABSTRACT 

In this paper we present the global baroclinic instability as a source for vigorous tur¬ 
bulence leading to angular momentum transport in Keplerian accretion disks. We show 
by analytical considerations and three-dimensional radiation hydro simulations that, in 
particular, protoplanetary disks have a negative radial entropy gradient, which makes 
them baroclinic. Two-dimensional numerical simulations show that a baroclinic flow is 
unstable and produces turbulence. These findings are tested for numerical effects by 
performing a simulation with a barotropic initial condition which shows that imposed 
turbulence rapidly decays. The turbulence in baroclinic disks transports angular mo¬ 
mentum outward and creates a radially inward bound accretion of matter. Potential 
energy is released and excess kinetic energy is dissipated. Finally the reheating of the 
gas supports the radial entropy gradient, forming a self consistent process. We mea¬ 
sure accretion rates in our 2D and 3D simulations of M = —10“® to — 10“^ Mq yr“^ 
and viscosity parameters of a = 10“^ — 10“^, which fit perfectly together and agree 
reasonably with observations. The turbulence creates pressure waves, Rossby waves, 
and vortices in the {R — f) plane of the disk. We demonstrate in a global simulation 
that these vortices tend to form out of little background noise and to be long-lasting 
features, which have already been suggested to lead to the formation of planets. 
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1. Introduction 

Protoplanetary disks appear to be a common feature around young stars (Beckwith & Sargent 
1993; Strom, Edwards, & Skrutskie 1993). They are thought to provide the material and the 
environment for the formation of planets (e.g. Lissauer 1993). Thus one needs to know the internal 
properties of such disks, such as the density, temperature and turbulence, in order to estimate the 
time scales of the formation process. These quantities are not directly accessible via observation, 
and so one needs a model for these disks to derive observable quantities like line emission and 
scattering efficiency for the light from the central star (Bell, Cassen, Klahr, & Henning 1997; 
Hartmann, Calvet, Gullbring, & D’Alessio 1998; Bell 1999). 

The basic idea of most models is that there is a process in these disks that transfers angular 
momentum radially outward, so that mass will flow radially inward (Liist 1952). Such a process 
might be turbulence (hydrodynamical or magneto-hydrodynamical) or self gravity (e.g. Larson 1989; 
Stone, Gammie, Balbus, & Hawley 2000). Independent of the source for the angular momentum 
transport, it can be parameterized by an effective viscosity v, which usually is scaled to the local 
sound speed c* and the pressure scale height Hd oi the disk by a dimensionless number called a 
(Shakura & Sunyaev 1973): 

V = aCsHo . ( 1 ) 

Using this simple description, it was possible to calculate the time evolution of disks (Lynden-Bell 
& Pringle 1974; Lin & Papaloizou 1980), the density and temperature structure, and the turbulent 
background as well as the laminar sub Keplerian mean component of the gas flow. The a models 
have influenced our understanding of the planetary formation process through their implications, for 
example, regarding the spatial distribution and collision rate of dust grains (Markiewicz, Mizuno, 
& Voelk 1991; Weidenschilling & Guzzi 1993; Dubrulle, Morfill, & Sterzik 1995; Klahr & Henning 
1997), the opening of gaps by giant planets (Bryden, Chen, Lin, Nelson, & Papaloizou 1999; Kley 
1999) or the radial drift rate of planets during their formation phas e (Goldreich & Tremaine 1980; 
Ward 1997). 

They are also used to explain the generation of FU-Orionis outbursts (Bell & Lin 1994; Kley 
& Lin 1999) as well as the measured accretion rates of T-Tauri stars (Hartmann et al. 1998). 

Despite the success of these models, a is still a free parameter for protoplanetary accretion 
disks. A barotropic Keplerian shear flow is in principle stable (Rayleigh stable) and does not 
develop turbulence per se despite the high Reynolds numbers. Dubrulle (1993), Richard & Zahn 
(1999), and Duschl, Strittmatter, & Biermann (2000) claim that the high Reynolds numbers will 
lead to turbulence in barotropic disks, but numerical investigations so far contradict this idea 
(Balbus, Hawley, & Stone 1996; Godon & Livio 1999a). Hydrodynamical turbulence has shown 
to transport angular momentum outward (Hawley, Balbus, & Winters 1999; Rudiger & Drecker 
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2001), but it is not able to sustain itself and usually rapidly decays in barotropic simulations. For a 
long time thermal convection in the vertical direction was thought to provide a turbulent transport 
of angular momentum in the outward direction (e.g. Cameron 1978; Lin & Papaloizou 1985). 
However analytical investigations of convection indicate inward transport of angular momentum 
(Ryu & Goodman 1992), a result that was later confirmed by numerical simulations (see below). 
For disks around black holes or in cataclysmic variables, magneto-hydrodynamical instabilities 
(e.g. magneto-rotational instability) seem to provide a reasonable amount of viscosity (Balbus 
& Hawley 1991). Relatively massive protostellar disks around young stars (> O.IM*) also have 
a proper mechanism to transport angular momentum via self gravity (Toomre 1964; Laughlin & 
Bodenheimer 1994). 

In the less massive protoplanetary accretion disks (~ O.OIM*) self gravity is only important 
at the outermost radii (Bell et al. 1997). At the same time the disk is so cool and contaminated by 
tiny dust grains, which capture most of the few free electrons, that magnetic fields have negligible 
influence on the bulk of the disk matter in the regions where planets are expected to form (Gammie 
1996). The sub-micron sized dust grains reduce the free electron density efficiently enough to 
prevent the possibility of the magneto-rotational instability, as the field lines can diffuse faster 
than the shear can tangle the field. Only in the vicinity of the star, where dust has evaporated 
and no planets can form, is there sufficient ionization for the development of magnetic turbulence, 
which can generate a k. 0.01 (Stone et al. 2000). In addition, if the optical depth above the 
disk is sufficiently small, cosmic rays and X-rays (Glassgold, Najita, & Igea 1997) could ionize the 
uppermost thin surface layers of the disk in the hypothetical absence of dust, which has led some 
authors to the speculation that between about 0.1 and 30 AU a disk consists of a non-accreting 
‘dead zone’ sandwiched between two active layers at the surface (Gammie 1996). 

Indeed, if there is no self gravity working, basically no ionization present and no shear insta¬ 
bility possible, the disk would develop no turbulence, transport no angular momentum, release no 
accretion energy and would basically cool down to the ambient temperature, while orbiting the 
star on a time-constant orbit. In such a disk, time scales for the growth of planetesimals would be 
much longer than in the standard scenarios, as the Brownian motion of the micron-sized dust is 
smaller and there is no turbulence acting on the dust in the cm- to m-size range. Sedimentation 
of the grains to the midplane could create a shear instability (Guzzi, Dobrovolskis, & Ghampney 
1993) which then generates turbulence. Some authors (Guzzi, Dobrovolskis, & Hogan 1996) argue 
that the solar nebula at one time must have been turbulent in order to explain the size segregation 
effects during the formation of chondrites. In any case we want to stress that it is of major relevance 
for planet formation to know whether protoplanetary disks are turbulent or not. 

In this paper we present numerical simulations of a purely hydrodynamical instability that 
works in accretion disks, namely a baroclinic instability, similar to the one responsible for turbulent 
patterns on planets, for example, Jupiter’s red spot and the weather patterns of cyclones and 
anti-cyclones on earth. Baroclinic instabilities arise in rotating fluids when surfaces of constant 
density are inclined with respect to the surfaces of constant pressure (e.g. Tritton & Davies 1985). 
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Vortensity, defined as vorticity per unit surface density, is not conserved as is the case in barotropic 
two-dimensional flows (Kelvin’s circulation theorem, see e.g. Pedlosky 1987), and vortices can be 
generated. 

Kippenhahn & Thomas (1982) studied the compatibility of thermal and hydrostatic equilibrium 
in thin radiative accretion disks but only checked axisymmetric flows. Cabot (1984) investigated the 
possibility and efficiency of a baroclinic instability for cataclysmic variables (CV) in a local fashion, 
concentrating on local vertical fluctuations. He found that this instability produces insufficient 
viscosity to explain CVs, but nevertheless enough viscosity (a w 10“^) for protoplanetary disks. 
Knobloch & Spruit (1986) also investigated the baroclinic instability due to the vertical stratification 
in the disk but found that it is not a reliable mechanism for angular momentum transport in thin 
disks. Ryu & Goodman (1992) considered linear growth of non-axisymmetric disturbances in 
convectively unstable disks. They used the shearing-sheet approximation in a uniform disk and 
found that the flux of angular momentum was inwards. Lin, Papaloizou, & Kley (1993) performed 
a linear stability analysis of non-axisymmetric convective instabilities in disks but allowed for some 
disk structure in the radial direction. Their disk model includes a small radial interval in which 
the entropy has a small local maximum of 7% above the background but with a steep drop with 
an average slope of iL ~ R~^ {K is the polytropic constant; see below). Such a situation could be 
baroclinically unstable, and in fact that region is found to be associated with outward transport 
of angular momentum. Lovelace, Li, Colgate, & Nelson (1999) and Li, Finn, Lovelace, & Colgate 
(2000) investigated the stability of a strong local entropy maximum (a factor 3 above the background 
with a slope K ~ in a thin Keplerian disk and found the situation to be unstable to the 

formation of Rossby waves, which transported angular momentum outward and ultimately formed 
vortices (Li, Colgate, Wendroff, Liska 2001). Sheehan, Davies, Cuzzi, &: Estberg (1999) studied the 
propagation and generation of Rossby waves in the protoplanetary nebula in great detail, but they 
had to assume some turbulence as a prerequisite. 

In contrast, here we investigate a baroclinic instability that arises from the general radial 
(global) stratification of the gas flow in accretion disks. Please keep in mind that global refers to 
baroclinic and does not necessarily imply that the instability has a global character, i.e. depends 
on the boundary conditions. Our motivation is based on the observation of positive Reynolds 
stresses in radiation-hydrodynamical three-dimensional (3D) simulations of thermal convection in 
protoplanetary accretion disks (Klahr & Bodenheimer 2000a). As convection is known to have the 
property of transporting angular momentum inward (Kley, Papaloizou, & Lin 1993; Cabot 1996, 
hereafter referred to as C96; Stone & Balbus 1996, hereafter referred to as SB96), we were surprised 
by this result and started an investigation to identify the special ingredient that would explain this 
contradictory result. We found that the size of the simulation domain, especially the azimuthal 
extent, influences the sign of the Reynolds stresses (Klahr & Bodenheimer 2000a,b). In order to 
accomplish a sufficient number of tests to isolate the crucial ingredient in our model and to show 
that artificial boundary effects are not contributing, we stripped the 3D radiation-hydro simulation 
down to a flat 2D (r — (p) disk calculation. The radially varying initial temperature (~ Entropy 
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~ K) in these calculations approximately reproduces the temperature distribution found in the 3D 
simulations. These tests show that indeed there is angular momentum transport outwards and that 
the lower-order azimuthal modes give the fastest growth rates in baroclinic simulations. 

In §2 we make an argument to show that protoplanetary accretion disks can be unstable to 
a non-axisymmetric global baroclinic instability as long as the source term of vorticity (baroclinic 
term) does not vanish. In §3 we explain the changes to the TRAMP code (Klahr, Henning, & Kley 
1999) that became necessary for the simulations presented here. Two results from 3D radiation 
hydro simulations are shown in §4. One model adopts an artificial heating term similar to that 
in C96 and SB96, while the other one is self-consistent in the sense that the only possible heating 
is a pure effect of compression and shock dissipation. In §5 we describe tests with the 2D flat 
approximation which prove numerically the existence of the instability and show the strength of 
the turbulence that develops as well as the Reynolds stresses that are produced. An initially 
isothermal (= globally barotropic) calculation which shows decaying turbulence demonstrates the 
reliability of our numerical investigation and proves that the instabilities observed in baroclinic 
models are not an effect of the shearing disk boundary conditions, or our numerical scheme. In §6 
we discuss the first global (2D) simulation of an accretion disk allowing for the baroclinic instability 
and find the interesting result that large vortices form. These vortices are long-lived high-pressure 
anti-cyclones with an over-density by a factor of up to four. It has been suggested that such vortices 
could be direct precursors of proto-planets, whose formation could be initiated by concentration of 
dust toward their centers (Barge & Sommeria 1995; Tanga, Babiano, Dubrulle, & Provenzale 1996; 
Godon & Divio 1999b). There is also the possibility that the over-dense regions could eventually 
undergo gravitational collapse (see Adams & Watkins 1995), but our present code does not allow 
for this process. In the last section (§7) we discuss our results. 


2. Stability considerations 

In this section we show that a baroclinic instability is plausible in a disk. The instability can 
only arise if there is an inclination between the density and pressure gradients (baroclinic term); 

V/9 X Vp / 0. (2) 

Here we analyze the situation using the standard polytropic equation 

P = (3) 

where p is the pressure and p is the mass density. We choose a value of 7 = 1.43 which is 
representative for a mixture of molecular hydrogen and neutral helium with typical abundances 
for protoplanetary accretion disks. If K is constant as a function of R, as used by various authors 
(e.g. Adams &: Watkins 1995, 1996; Goldreich, Goodman, & Narayan 1986; Nautha & Toth 1998; 
Godon & Divio 1999a,b), the baroclinic term vanishes and vorticity is conserved in plane parallel 
flows. 
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In contrast to previous numerical work we choose a radially varying entropy K{R) (§5.2) for 
the initial state of the 2D simulations. This is the best way to mimic the density and temperature 
distribution that arises in the 3D radiation-hydro calculation. In other words any realistic disk has 
a radial entropy (5) gradient {p/^ const). This can be seen from a simple argument, where we 
only assume that the aspect ratio of the disk {H^/R) is constant for a range of radii. The pressure 
is never only a function of local density but is also a function of local gravity, both of which change 
with radius; 

p ~ ~ pn‘^R^ ~ pR-^ ~ (4) 

This result holds for any accretion disk, as long as H^/R is constant. The exact dependence of 
K on R nevertheless depends on the specifics of a given simulation. For example in a thermally 
convective region a density distribution of p oc R~^ and a temperature distribution of T oc R~^ are 
typical. Thus it follows 

p oc Tp oc K{R)p"' => K{R) oc ( 5 ) 

In a more general fashion, if p oc R~^p, T oc R~^'^, and K oc R~^^ this equation reads; 

f^K = f3T + Pp (1 - 7 ) • (6) 

It obvious, that there are certain stable profiles. For if Pt = Pp{l — 1) then Pk = 0 (e.g. for 
Pp = 1 => Pt = 0.43). But our 3D radiation-hydro calculations do not show this particular profile 
(see § 4). Thus Pk = 0.57 for the initial state in our baroclinic simulations. 

It follows from equations (2) and (6) that 

Vp X V{K{R)p^) ~ -Pk {d^p) , (7) 

which is always non-zero if there is an azimuthal density fluctuation and a non-zero Pk- On the 
other hand, isothermal disks are always barotropic by definition. 

Baroclinic instabilities have been widely studied in the context of meteorology and oceanog¬ 
raphy (e.g. Tritton & Davies 1985, Pedlosky 1987). There are some theoretical models as well as 
some laboratory experiments which can help to understand the basic mechanism of this instability 
(Sakai, lizawa, & Aramaki 1997). The onset of a baroclinic instability lies in the radial entropy 
gradient. This entropy gradient can result from the temperature gradient between the equator and 
north pole of the earth. It can also be realized in an experiment in which a fluid is contained 
between two concentric cylinders - one is heated and the other one is cooled down. In the absence 
of rotation thermal convection will occur in the radial direction. This would lead to one large scale 
convection cell between equator and north pole (Hadley cell, Hadley 1735) as well as to thermal 
convection between the cylinders. But the rotation of the earth (or the rotation of the cylinders) 
inhibits this convection as it allows only for concentric flows and no convective heat exchange occurs 
in the radial direction. But if the entropy gradient is strong enough, the concentric flow becomes 
unstable and begins to meander. This is the baroclinic instability. As the water meanders, it is 
once again able to transport heat radially. The same happens in the earth’s atmosphere, where the 
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baroclinic instability leads to meteorological effects at mid-latitudes like the meandering jet stream 
and the high and low pressure systems (anti-cyclones and cyclones), determining our daily weather. 

Protoplanetary accretion disks are believed to have thermal convection perpendicular to their 
midplane (Cameron 1978; Lin & Papaloizou 1985). The disks do not fulfill the Schwarzschild crite¬ 
rion for vertical stability at a wide range of radii. Also for the radial stratification a Schwarzschild 
criterion can be formulated which indicates that if rotational effects were negligible there would be 
radial thermal convection occurring (e.g. convective zones in stars). 


The Schwarzschild criterion for stability of the radial stratification is actually never fulfilled in 
our baroclinic disks, as the radial component of the adiabatic temperature gradient (see Kippenhahn 
& Weigert 1990) is 

PdT I 

V... = -^ = l-- = 0,3 (8) 

while the absolute temperature gradient in the radial direction is 


PdT _ Pt 
T dP /3t + Idp 


= 0.5. 


(9) 


Thus Vad < V. Furthermore the Brunt-Vaisala frequency in the radial direction can be calculated 
via (see Kippenhahn & Weigert 1990): 


N\r) = -^{Vad-V). (10) 

Pd 

Here g is the radial component of gravity (not its absolute value as in Kippenhahn & Weigert). 
At a local radius Rq centrifugal force and gravity may cancel each other. But if a parcel of gas is 
radially perturbed in a viscosity free disk, then it conserves its given angular momentum Hoi?o- 
will therefore feel an effective gravity of 


g = nl {R^R-^ - RlR-^) , 

which can be expanded for small deviations a = R — Rq into 

g = —Q^a. 


( 11 ) 


( 12 ) 


Radial gravity follows locally basically the same dependence as the vertical component of gravity. 
Thus it follows 

Ar2(r) = -0.211^ (13) 

which is obviously negative and hereby convectively unstable. Interestingly, a positive gradient in 
entropy also leads to this situation, but in this case the flow would be unstable to perturbations 
with a < 0. 


Nevertheless, the rotation prevents the disk from being convectively unstable in the radial 
direction, which can be seen from the baroclinic version of the Rayleigh criterion (e.g. Solberg- 
Hpiland criterion) 


1 dR^n^ 
W dR 


dlnP 

dR 


(V - Vad) > 0 


(14) 



which in the case of a Keplerian protoplanetary accretion disk is: 


> 0. (15) 

With p ~ pT ~ R~^ it follows 

1 - 0 . 4 :^) >0, (16) 

which gives a value of about 0.9611^, which is a sufficient condition for stability if only axisymmetric 
perturbations are permitted. If non-axisymmetric perturbations are allowed it is only a necessary 
condition for stability. Such a situation is a typical onset to the formation of non-axisymmetric 
baroclinic instabilities (meanders, vortices) as we discussed before (see Tritton &: Davies 1985). 

A detailed linear or non-linear, possibly local, stability analysis of the global baroclinic insta¬ 
bility in the context of accretion disks still has to be performed to get a sufficient criterion on the 
instability. Such an analysis should not be confused with the global stability analysis of a local 
baroclinic instability as in Lovelace, Li, Colgate & Nelson (1999) and Li, Finn, Lovelace, & Colgate 
( 2000 ). 



3. Changes in TRAMP 

The code Tramp (Three-dimensional RAdiation-hydrodynamical Modelling Project) was 
introduced and extensively discussed by Klahr et al. (1999). The equations used here are the same 
as in that paper. However, two additional features have been added to the code, which are crucial 
for the work presented in this paper. 


3.1. Shearing Disk 

As already mentioned in Klahr et al. (1999), rigid wall boundaries in the radial direction 
can affect the flow-pattern and are not the perfect choice for simulations of disks. SB96 use the 
ZEUS code described in Stone & Norman (1992) in the shearing sheet approximation (Hawley 
et al. 1995). This approximation is only possible for pseudo-Cartesian coordinate systems which 
themselves strongly limit the scale of the computational domain, i.e. the box size has to be small in 
comparison to the local radius. A direct consequence of the pseudo-Cartesian coordinates is that 
there can be no global radial gradient of the density, temperature and entropy. Thus, local studies 
of the global baroclinic instability are impossible. 

In order to overcome this problem. Tramp calculates on a spherical (i?, 6,4>) grid; thus large 
portions of the disk can be simulated, and one is only restricted by the computational effort required. 
We put forward a new set of radial boundary conditions that we call shearing disA: boundaries, which 
in fact correspond to a shearing sheet for spherical coordinates and disks. In the shearing sheet 
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approximation one uses simple periodic boundary conditions, i.e. the values for the ghost-cells are 
determined from the value of cells on the opposite side of the computational domain in a simple 
fashion [e.g. for the density /3(0) = p{NX — 1), where NX is the maximum number of zones in 
a given direction]. In the shearing disk approximation one makes two assumptions: first, that the 
mean values of each quantity follow a power law scaling with the radius R, which is usually the 
case in steady-state accretion disks for a certain range of radii, and second, that the fluctuations 
are proportional to the mean value. Thus it follows for the two inner and two outer ghost cells: 


pi0,J,K) 

= p*{NX-l,J,K){R{C)/R{NX-l))-!^”, 

(17) 

p{l,J,K) 

= p*{NX,J,K){R{l)/R{NX))-^^, 

(18) 

p{NX + l,J,K) 

= p*(2,J,A)(4^(iVA + l)/i^(2))-^^ 

(19) 

p{NX + 2, J, K) 

= p*{‘i,J,K){R{NX+ 2)/R{‘i))-^p. 

(20) 


The {*) indicates the implementation of the global shear in a manner similar to that given in SB96 
and C96. With the assumption that the inner and outer boundaries each move with the local 
Keplerian speed, it follows that there is a mean angular offset between inner and outer grid at time 
t oi 64> = “ ^R(NX) ~ ^RiNX+i))^/"^- This offset translates into an integer offset 

dK = ABS{5(j)/d(j)) defined via the angular width of a grid cell d(j) and two interpolation factors 
Cl = mod{5(j),d(j))/d4> and C 2 = 1 — Ci: 

p*{K) = C2p{K + dK) + Cip{K + dK- 1). (21) 

The value for (dp is 1 in all 3D simulations as we are performing calculations around R = bAU 
where the surface density is almost constant (/3s = 0). For the azimuthal frequency {g = (p = to) 
= 1.5 is obvious, = 1.0 in the 3D simulations for the temperature follows from the assumption 
of a constant relative pressure scale height Hp/R = const, which is necessary to make inner and 
outer vertical stratification fit in polar coordinates. The 2D simulations use (St = 0.0 in the 
isentropic case (model 2) but (3t = 1-0 in the baroclinic unstable case (models 3,4 and 5). Model 6 
uses open boundaries rather than shearing-disk conditions. Finally fdg = 1.5 for the polar velocity 
component {g = 9) assumes the fluctuations to be proportional to the local speed of sound. For 
the radial velocity Ur, Pu = — f3p is chosen to conserve the mass flux over the radial boundaries 

in a spherical coordinate system where pUrR^ = const. It was tested that kinetic energy flux 
and momentum flux are not artificially amplified by the boundary condition (see model 2). We 
damp the radial velocity component of the inner and outer four grid cells by a factor of fd = 0.05 
each time step (ur := Ur{1.0 — fd)) in order to prevent the model from creating artificial resonant 
oscillations in the radial direction. These radial waves can even occur in ID radial models if there 
is no damping. Model 2 (see below) is our basic test that the treatment of our boundary conditions 
does not lead to a numerical instability. Angular momentum is not conserved, which allows the 
possibility of a net transport of angular momentum into or out of the computational domain by 
turbulence, and thus of driving an accretion process or suppressing it, depending on the direction 
of the transport. 
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This new method was extensively tested and will be the subject of an upcoming paper where 
its properties will be compared to other choices of boundary conditions under various physical disk 
situations. 

Using this new numerical method we combine the virtues of local simulations (good resolu¬ 
tion and moderately large time-steps) with the possibility of treating global systematic gradients, 
previously exclusively reserved for global simulations. 


3.2. Reynolds Averaging 


As an addition to the earlier Tramp version there is now the possibility to measure the trans¬ 
port properties of the turbulence, i.e. the correlations of the fluctuations especially for the radial 
transport of angular momentum. A similar method was already used in C96 as well as in SB96. 
This turbulent stress tensor T is then scaled by the square of the sound-speed, representing the 
pressure at the midplane. The component corresponds to the angular momentum transport 
which corresponds to the classical a value, or in other words, an a is a measure of how much vis¬ 
cosity would be needed in order to simulate the turbulent angular momentum transport by viscous 
diffusion of angular momentum in a laminar flow. As we also investigate other components of the 
stress tensor we define: 

< u'ApUr)' > 

(X ^r<h —9 

pci 



with 


< u'^{pUr)' > = < pU^Ur > —U^ < pUr >, 


(23) 


where <> or a bar represents a spatial and time average. We use the symbol Aij because the 
symbol Onf, is generally used for the helicity of the velocity field. 


A detailed discussion of the measurements of the Reynolds stresses including Are and for 
compressible turbulence shall be discussed in an upcoming paper (Klahr & Rudiger in preparation). 
In a similar manner we measure the strength and isotropy of the turbulence via the determination 
of the on-diagonal terms of the stress tensor Trr, Tge and which, scaled by the square of the 
sound speed, are three different indicators of the square of the turbulent Mach number: 


and 


^ _ <u'r^> 


^ee — - 5 — 


< > 



In the same way we determine the mean density fluctuation in the flow by 
M = Arr + Agg + A^^ to be the Mach number in our models. 


(24) 

(25) 

(26) 

^ . We define 
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4. 3D Radiation-Hydrodynamical Simulations 

The inviscid 3D simulations of thermal convection in disks (Klahr & Bodenheimer 2000) showed 
radially outward directed angular momentum transport if the section of the disk was large enough. 
In the following we show that the angular momentum transport is not dominated by the vertical 
thermal convection but by the hydrodynamical turbulence in the R — (p plane of the disk which 
results from the radial entropy gradient self-consistently introduced by solving the energy equation 
in our radiation hydrodynamical code. 

As a continuation of our previous work (Klahr et al. 1999), we have performed simulations of 
3D chunks of protoplanetary accretion disks. We have tested the influence of viscosity, artificial 
heating and boundary conditions extensively. As highly viscous flows tend to become axisymmetric 
(C96; Klahr et al. 1999), we were especially interested in simulations with high Reynolds numbers 
as in SB96. In model 1 we switch off the application of the viscous forces in the momentum equation 
but keep the heating, corresponding to a = 0.01, from the local viscous dissipation function (<1>) in 
the energy equation 

= (VT) u, (27) 

where T is the viscous stress tensor. For the viscosity v we adopt v = Kepler- For more 

details we refer to Klahr et al. (1999). In an additional simulation (model lb) we forgo the artificial 
heating and show that a self-consistent heating of the disk via the dissipation of compression waves, 
especially shocks, and the release of gravitational energy can lead to the same baroclinic effects as 
the artificial heating in model 1. Artificial heating simplifies the numerical effort, as it stabilises 
the vertical and radial structure of the disk and thus allows for longer integration times for the 
mean values (see Table 1). 

Artificial von Neumann & Richtmyer (1950) viscosity (see also Stone & Norman 1992) for 
the proper treatment shocks is included, as is flux-limited radiative transport with a simple dust 
opacity (k = 2 x 10“^T^ cm^ g~^)- The goal was to investigate the influence of our more global 
approach, especially a wider azimuthal range, on the simulation results. 

The particular models we present here (model 1 & IB; see Table 1) cover a radial range from 
3.5 AU to 6.5 AU, a vertical opening angle of ±7°, and an azimuthal extent of 90°. They use only 
20 grid cells in the vertical direction but 51 in radius and 60 in azimuth. A follow-up paper will 
deal entirely with 3D models at various resolutions. At the moment it is sufficient to say that the 
results of model 1 are typical and remain the same at much higher resolution. The benefit of model 
1 is that it can easily be evolved for 68 orbits at the outer edge of the simulation within sixteen 
CPU hours on a Cray T90. In future papers we will discuss the influence of artificial heating in 
contrast to self-consistent heating due to shock dissipation in greater detail. Calculations as model 
IB indicate at least that numerical dissipation might be a problem as it removes kinetic energy 
from the system without heating the gas. Thus there is less heating to maintain the initial vertical 
and radial temperature structure of the disk. Higher resolution models will have to investigate the 
effect of numerical dissipation. 
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Figure 1 shows the three-dimensional rendering of the vertical mass flux distribution after 
model 1 approached a quasi-steady state. The azimuthally extended convection cells are slightly 
twisted and not completely axisymmetric. In Figure 2 the normalised Reynolds stresses and corre¬ 
sponding Mach numbers of the turbulence for model 1 are given, time-averaged over 68 outer orbits. 
The radial mean value for the angular momentum transport corresponds to a = Ayip = 2 x 10“^, 
a quite high value, and the tendency towards positive values is obvious. The radially varying 
Reynolds stresses reflect their connection to locally and temporarily confined flow structures and 
events. Only a much longer time for averaging, as is done in the 2D simulations, can produce a more 
homogeneous distribution. Based on the measured stresses one can estimate a viscous time-scale 
of ~ a~^ = 500 orbits, while for technical reasons only 68 are possible in one run. The on-diagonal 
terms of A, which characterise the strength and isotropy of the turbulence, show that the fluctua¬ 
tions in radial and azimuthal velocity are almost sonic at some radii, and that the fluctuations in the 
vertical (polar) velocities are more than an order of magnitude smaller. This difference reflects the 
fact that (1) the radial and azimuthal velocity fluctuations are not driven by the vertical thermal 
convection, and (2) there is no isotropic (3-D) turbulence present. The resulting Mach number of 
the turbulence, the combination of the strengths in three dimensions, is close to one. The relative 
density fluctuations are in the mean as strong as 10%. 

In contrast to this result, with positive values of the Reynolds stress, our corresponding axisym¬ 
metric 2D simulations, with the same physical assumptions and parameters as in model 1, always 
showed clearly negative Reynolds stresses. The dissipation-free axisymmetric Euler equation would 
predict precisely zero average stresses, which is a consequence of strict angular momentum con¬ 
servation, but radiation hydro calculations are by nature not dissipation free even if one has no 
viscous forces. Thus the negative a values are a result of the heating and cooling of the gas. 

Looking at the flow-pattern and density distribution in the R — (/>-plane (Fig. 3), we do not 
see the convection cells from Figure 1. What we see are vortices, vorticity waves (e.g. Rossby 
waves, baroclinic waves), and pressure waves propagating in the flow. These are the sources of 
the positive Reynolds stress. They are counteracting the negative Reynolds stresses generated by 
the thermal convection. In the turbulence terminology one would possibly identify geostrophic 
turbulence (irregular waves; see e.g. Tritton & Davies 1985). This can be seen in the spectral 
density distribution of the flow in the azimuthal direction (Fig. 4). Only at the smallest resolved 
scales (i.e. meso scale) the slope approaches the Kolmogorov spectrum At larger scales 

the slope is even steeper than for earth atmospheric geostrophic turbulence {k~^) which indicates 
how strong small-scale turbulence is being pumped into the smallest possible wave numbers. In 
this calculation the minimal k is four as we calculate only a quarter of the disk. However 360° 
full circle simulations (see model 6) indicate that A: = I is the mode where the energy will pile 
up. The dissipation in this calculation thus does not only occur at the smallest scales as in 3D 
incompressible turbulence but at all scales, especially at the large scales where shocks form. The 
idea to drive the accretion process via shocks was already suggested by Spruit (1987), but in his 
barotropic models the shocks could not form without the presence of an external perturber. 
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The Mach number of the flow is the smallest in the midplane (compare M^ax = 0.6 in Fig. 3 
with M = 0.5 in Fig. 2). Velocities rise, eventually beyond the sound speed, with height above the 
midplane. 

From this result several questions arise. (A); Why did neither Stone & Balbus (1996) nor 
Balbus, Hawley & Stone (1996) observe these positive Reynolds stresses and violent turbulence? 
(B): What is the source for the vorticity and turbulence in the R — 4> direction? (C): Especially, 
how can we be sure that we are not observing boundary effects? We will explicitly answer these 
questions in §7. 

The strategy in order to identify the source for turbulence and angular momentum transport 
was to simplify our model step by step by removing physics from the simulation and to wait for 
the turbulence and Reynolds stresses to disappear. The first effect we found was that whenever we 
decreased the azimuthal extent of the computational domain below a critical value of ~ 15 degrees, 
the Reynolds stresses switched to negative values, as in 2D axisymmetric calculations (Klahr & 
Bodenheimer 2000b). As the SB96 calculation only covers roughly 1 degree, this was a first hint 
that our findings might result from the more global simulation. 


4.1. No artificial heating 

The artificial heating in model 1 has the benefit that there is a well defined laminar equilibrium 
state for the disk. This is ideal for any stability investigation as well as for the setup of a numerical 
model. 

But one can argue that this artificial heating is inconsistent with the lack of dissipation of 
kinetic energy. Consequently, the findings with model 1 could result purely from the artificial 
heating, or, even worse, from the specific shape of the dissipation function (see Eq. 27). Thus, 
we repeated the simulation of model 1, but switched off the heating completely. This simulation 
is numerically more complicated and more expensive than model 1 as the initial state is far from 
equilibrium and the disk undergoes strong fluctuations in density and temperature. Only 14 orbits 
could be performed during a 10-hour run on the Cray T90. The disk vertically shrinks by 20 percent 
before it reaches a new self-consistent state with a midplane temperature of 40 K, down from 150 
K in the initial state. 

If the temperature changes by a factor of almost 4 (150 K ^ 40K) one could expect the disk 
height (~ pressure scale height ) to change by about a factor of 2. But the factor of 2 in scale 
height is only true for vertically isothermal disks. Our disks have a vertical temperature profile 
that becomes flatter as the temperatures decrease because the opacities also decrease (proportional 
to K oc T^). This means that the high-temperature disk has temperature, pressure and density 
profiles that fall off much more steeply than the pressure scale height estimated for the midplane 
temperature might suggest. The cooler disk on the other hand is closer to an isothermal disk and 
thus does not shrink in proportion to the root of the midplane temperature. 
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The vertical grid structure was readjusted during this cooling: the initial vertical opening angle 
was ±7°, but during the vertical shrinking the gradients above the disk became too steep for the 
numerical technique to work. Thus we constructed a new vertical grid ±5.5° also with 20 grid 
points and interpolated the values linearly from the old onto the new grid. 

The newly gained state of 40 K was then maintained for more than 100 orbits (see Fig. 5, which 
is similar to Fig. 3). The disk turbulence was transonic (see Fig. 6) and the heating was highly 
non-homogeneous in space and time, in contrast to model 1 where it was assumed to be smoothly 
distributed. But nevertheless the general result was the same: the disk develops a radial entropy 
gradient, turbulence and vorticity. In addition, the accretion of matter which results from the 
positive Reynolds stresses feeds energy into the system which is radiated away after the turbulence 
is dissipated in shocks. 

These self-consistent simulations will have to be continued in our future work. The simulations 
are numerically much more unstable than the heated or polytropic models. But in the end only 
they can provide realistic predictions for density, temperature, turbulence, and the evolution of 
accretion disks. 


4.2. Radial entropy gradient 

Another glance at Figure 3 shows that temperature (contour-lines) and density (colors) are not 
perfectly aligned. Thus also pressure and density gradients do not point in the same direction, a 
clear indication of a baroclinic flow, where the entropy decreases with radius. To test the relevance 
of this possible instability for the generation of the observed turbulence we removed the vertical 
structure of the disk and got rid of the radiation transport. Thus the 3D density p was replaced 
by the surface density S and the 3D pressure p by the vertically integrated pressure P = ST (in 
dimensionless units). In Figure 7 we plot the K ~ Tp^~'^ values in the midplane (before model 1 
became turbulent) as K/Kmax- If one measures this K value in the turbulent state of the disk, 
the fluctuations of K are too strong, which hides the general trend. The next section describes the 
results of using such an entropy distribution in a two-dimensional (R, 4>) disk. 


5. 2D Simulations 

The flow characteristics in the 3D simulations, e.g. the turbulence cascade and the small vertical 
velocity fluctuation, indicated that the disk instability must be a 2D effect. In order to verify this 
thesis we removed the vertical structure and performed 2D (i? — (f>) simulations. 

Both of the following models were first developed as ID radial axisymmetric models with 
density slope p ~ R~^ {p[lAU] = 10“^® g cm“^), which corresponds, in the 2-D flat disk, to S = 
const (~ 300 g cm“^). The initial temperature was then chosen to be either constant (T = Tq) as 
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a function of radius or to follow a power law (T = Tq {R/Ro)~^). The value for Tq was adjusted 
to create a local pressure scale height of H/R = 0.055 at the mid-radius Rq = 5AU, to match that 
in the 3D radiation hydro calculations. For obvious reasons this pressure scale height condition is 
only fulfilled at Rq for the constant temperature case but for all radii in the T oc R~^ case. In 
the constant temperature model H/R varies as which means that the relative pressure scale 

height increases slightly with radius. Both models use an identical computational domain with 
radii between 4.0 and 6.0 AU and an azimuthal extent of 30°, with a moderate resolution at 64^ 
grid cells. The boundary conditions are identical to those of model 1: a shearing disk plus damping 
of the radial velocities in 4 grid cells near the inner and outer boundaries. The artificial viscosity 
is the same as in model 1. Both models are initially kicked by a random density perturbation with 
1% amplitude; p := p * (0.99 -|- 0.02 * RAND), where RAND is a uniformly distributed random 
number between 0 and 1. Afterwards the system was allowed to evolve freely. We advect the 
internal energy and calculate the PdV work arising from compression and shock dissipation in the 
same manner as in the 3D calculations. For the time being there is no radiative cooling, because 
cooling would immediately change the temperature profile and we are interested in the behavior of 
given temperature distributions. Models with radiative cooling shall be dealt with in future work. 


5.1. Model 2: radially constant entropy: (T = const) 

This model developed no turbulence. The initially induced turbulence decayed rapidly with 
an e-folding time for the kinetic energy in i? of ~ 4.3 orbits (see Figs. 8 and 9). Nevertheless, the 
initially introduced turbulence always transported angular momentum outward, never inward (see 
Model 2B). Figure 8 shows the surface density distribution and flow field after about 90 orbits. 
The surface density is almost constant (299.955 < S < 300.125 g cm“^). The arrows indicate that 
there is still small noise in the velocity field with a Mach number of about 10“^. In Figure 9 the 
evolution of the total kinetic energy in both coordinates, and the enstrophy of the flow are plotted, 
in units of the initial values. The overall mean mass accretion rate is given in solar masses per 
year. Here and in the following models, a negative value of M indicates accretion towards the star. 
Enstrophy is the integral square of the local vorticity. One clearly sees that the kinetic energy 
as well as enstrophy are decaying as is expected for a barotropic disk. Kinetic energy in cj) and 
enstrophy are calculated from the deviations from the Keplerian profile. The residual finite values 
of these quantities result from the laminar (not strictly Keplerian) steady state flow. The resulting 
mean turbulence, measured by the strength of the components of the stress tensor (see Fig. 10) is 
also very low. 

We also tested the effect of stronger initial perturbations. For model 2B we gave model 2 an 
initial random density fluctuation of ±50%. Still the turbulence rapidly decayed as expected. But 
we have to stress that the turbulence in this barotropic simulation had already the property of 
effectively transporting angular momentum outward. The Reynolds stresses during the first orbit 
(see Fig. 11) are positive {< a = A^^ >~ 3.0 x 10“^), nevertheless rapidly decreasing afterwards. 
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We conclude that positive Reynolds stresses are a common feature of hydrodynamical turbu¬ 
lence in disks, and not exclusively a result of baroclinic simulations. This is also a result of Hawley 
et al. (1999), which shows that positive Reynolds stresses are correlated to the decay of turbulence. 
Only thermal convection, which must not be confused with isotropic turbulence, has a tendency of 
transporting angular momentum inward. Nevertheless, in barotropic disks hydrodynamical turbu¬ 
lence can not sustain itself and decays. 


5.2. Model 3: radially varying entropy: (T ~ R 

Using a different initial state leads to a completely different result. The flow becomes turbulent 
within a few orbits, with kinetic energy rising exponentially until saturation due to shock dissipation 
occurs. Model 3 (Fig. 12) shows some similarities to model 1 (Fig. 3) even though we are comparing 
a 3D and a 2D simulation. After 100 orbits the surface density shows deviations from axisymmetry 
(286 < S < 318 g cm“^). The velocities are 10^ times stronger than the ones in model 2. 

The initial instability grows quickly (Fig. 13). The kinetic energy in the radial direction grows 
by a factor of 10^ within 40 orbits. This translates into a characteristic growth time of ~ 5.0 orbits. 
The other components, especially the enstrophy, grow more slowly but do not saturate as quickly 
as the radial kinetic energy (which is also slightly damped for numerical reasons [see §3.1]). The 
strong rise of enstrophy (~ vorticity generation) is an indication that a baroclinic instability is at 
work. 

The angular momentum transport and also the on-diagonal components of the stress tensor 
(Fig. 14) are orders of magnitude stronger than in model 2. They are weaker than in models 
1 and IB, which may result from the smaller computational domain, especially in the azimuthal 
direction, artificially limiting the wave numbers of the instability. The angular momentum transport 
(a = A,.^) is in the mean as strong as 1.5 x 10“^. This value should be contrasted to the strength of 
the turbulence itself. An A^r and A^^ of ~ 10“^ correspond to a turbulent Mach number as strong 
as 0.05. In isotropic turbulence this value could be used for a mixing length model and would 
then predict an alpha-viscosity of a = 2.5 x 10“^ — 5 x 10“^, depending on the estimates for the 
“typical eddy size”. The lower value of a that we obtain indicates that even when strong turbulence 
develops it does not automatically generate strong Reynolds stresses, which again is evidence for 
non-isotropic turbulence. The energy of this turbulence is drawn from the entropy background 
which can maintain its profile due to the accretion of mass. The measured mass accretion rate of 
M ~ —2.0 X 10“® Mq yr“^ is in rough agreement with the expectation from a viscous accretion 
disk model with a surface density of ~ 300 g cm“^ and an a parameter of 1.5 x 10“^. 

The turbulence saturates due to dissipation from shocks and stays at a high, almost sonic, 
level. This level is on the long term quite variable. 

Our simulation demonstrates that non-magnetic turbulence can drive outward angular mo¬ 
mentum transport and is maintained itself by the resulting accretion process. Such a conclusion 
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has been in doubt for a long time, as no instability mechanism has been convincingly demonstrated 
to work in Rayleigh-stable and non-ionized disks. 

We checked our findings on model 3 by redoing the simulation at twice the resolution (model 
4; 128 X 128 grid cells) and obtained general agreement (see Fig. 15) with a mean a = Ar^ slightly 
larger than that of model 3. A detailed study of the influence of resolution on our results is a 
current subject of our investigations and will be part of a future paper. 

For model 5 we used also a resolution of 128 x 128 grid cells, but this time we used a compu¬ 
tational domain which was essentially twice as wide in both directions. It spans radii from 3 AU 
to 7 AU and 60° in the (/>—direction. The numerical resolution is thus about the same as in model 
3. In Fig. 16 we plot the Reynolds stresses, which are stronger than in the previous models with 
a = Artj, ~ 6 X 10“^. This larger scale simulation showed also the formation of a vortex (see Fig. 
17). It is an anti-cyclonic prograde vortex with higher density and pressure than the back ground. 
The vortex emits spiral waves in the ambient accretion disk. 

As model 5 is so far the best compromise between resolution and angular width, it is well 
suited for the examination of the spectral density distribution (see Fig. 18) and comparison with 
model 1 (see Fig. 4). Here one clearly sees the transition between 2D and “3D” characteristics 
(even though the calculations are actually 2-D) of the spectrum at a wave number of about 60, 
a value which was so far observed in all simulations. This transitional wavenumber is known as 
“Rhines blocking wavenumber”; roughly speaking, for m = /c > 60 energy only cascades down 
to smaller scales while for /c < 60 it inversely cascades up to larger scales. The wave number 
k = 60 corresponds to about 6 degrees in azimuth or a tenth of the radius, which is almost two 
pressure scale heights. Structures that are smaller than this are not influenced by the rotation and 
apparently behave like 3D turbulence i.e. energy cascades down to smaller scales to be 

ultimately dissipated. Structures bigger than this feel the Coriolis force and thus the rotation and 
shear. Consequently they are forced to have 2D characteristics and an inverse energy cascade with 
k~^. The same behavior can be observed in the earth’s atmosphere, where small structures up to 
500 km in extent (i.e. local winds) behave according while larger structures (like hurricanes) 

follow a k~^ distribution. 

The studies on the influence of the computational domain will have to be continued in future 
investigations. The same is true for parameter studies on the initial distribution of density (as 
well as 7 and temperature) in the disk models and the influence on a. The fact that we found a 
working hydrodynamical instability for disks which generates Reynolds stresses with the right sign 
and reasonable a-values, seems to support the paradigm of an viscously driven accretion process. 
Considering the limitations of our numerical method, we could argue that a better method or higher 
resolution could also lead to even stronger turbulent viscosity, as only the sound speed seems to be 
a natural upper limit for the velocity fluctuations. 

Other flow characteristics besides the Reynolds stress are very strongly developed. First, the 
2D-turbulence itself produces a significant turbulent (r.m.s.) pressure which is not accounted for in 
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a—models. Second, the vorticity and the deviation from the mean Keplerian profile are significant 
features, which also cannot be handled in diffusion models. And finally even though a— models 
might assume that there is some underlying non-axisymmetric turbulence, they never would predict 
the formation of large scale flow features like long-lived and growing vortices of several scale heights 
diameter. Even if such features were assumed to be present initially, a— models would smear them 
out and not amplify them (Godon & Livio 1999b). 

We have now demonstrated the capability of the global baroclinic instability to form non- 
axisymmetric structures and vortices in the local simulations with periodic boundaries. The next 
question to be investigated concerns what would happen if the global shear were not fixed by the 
radial boundary conditions but were free to evolve. Thus we set up a global baroclinic simulation 
to see how a real disk would behave and to what extent it could still be described by an a—model. 


6. A global 2D simulation 

Our global simulation shows that all the findings from the local models (model 3,4 and 5) can 
be confirmed. Actually the Reynolds-stresses are even stronger than in the local simulations, as 
smaller wave numbers in the azimuthal direction can be resolved. The wave number m = 1 seems 
to be the preferential mode of the instability. A further discovery was the self-consistent formation 
of long-lived anti-cyclonic vortices as a direct result of the global baroclinic instability. This may 
have major relevance for the formation process of planets. 

The simulation (model 6) covers the entire 360° of the circumference of the disk and a radial 
section between 1 AU and 10 AU. The grid measures 128 x 128 grid-cells, which are radially 
logarithmically distributed. We can thus resolve azimuthal wave numbers between 1 and 64. The 
boundary conditions in the radial direction were changed from periodic to simple non-reflecting 
outflow conditions (vanishing gradients), not allowing for inflow. As a result of this change we 
also could drop the damping of the radial component of the velocity close to radial boundaries. 
The density distribution again was p oc R~^ (constant S 300 g cm“^), and the temperature 
distribution was T oc thus we have a baroclinically unstable situation as in model 3, which 
results from H/R = 0.055. The model was first run into a stable ID axisymmetric state, where 
the residual velocities were less than 10“^ cm s“^. Without a symmetry-breaking instability and 
turbulence generation, this disk can not evolve and would stay perfectly laminar forever, as in the 
“dead zone” described by Gammie (1996). 

The initial density distribution was then perturbed by random noise of amplitude only 0.1%. 
The initial state is practically axisymmetric. Figure 19 illustrates the evolution of the flow in 
two space dimensions, over the full 360°. After the first 1 orbit (30 yr at 10 AU; Fig. 19a) only 
little structure has evolved. But with time a prominent anticyclonic vortex forms, which reflects 
the assumption that m = 1 is the preferred mode. Intermittently a second vortex also forms, and 
we assume that their number is just limited by the narrowness of our disk and a lack of matter. 
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The vortex grows in mass and propagates radially outward, possibly as a result of the gradient of 
background vorticity and the fact that anticyclonic vortices are a local sink for angular momentum 
in the global vorticity field. Anyway, as there was never a radial drift of vortices in barotropic 
flows reported, this effect might also be linked to the baroclinic features of the flow. A detailed 
investigation of this effect still has to be performed. 

Figure 20 shows the Reynolds stress and turbulent Mach number averaged over the orbits 
from 430 to 500. We see that the Reynolds stress does not disappear as we remove the “shearing 
disk” boundaries. The stresses are comparable in the main part of the disk to those in models 3, 
4, and 5 and only deviate at the physical edge of the disk where density and sound speed drop 
by orders of magnitude. We can conclude that the “shearing disk” boundary condition does not 
significantly affect the results. 

Figure 21 shows the situation after about 10^ years in real Cartesian coordinates to give an 
impression of the global nature of the simulation. Figure 22 shows the flow pattern in more detail 
in (i?, (p) coordinates at the end of the simulation. A huge vortex has formed which has a factor 4 
over-density with respect to the ambient disk and a factor 2 over-density with respect to the initial 
local surface density. It is a high-pressure anticyclone that has the property of collecting solid 
material in its center (Tanga et al. 1996; Godon & Livio 1999b). At the same time the over-dense 
blob inherits a substantial fraction of the disk gas, which is not confined by self gravity but only 
by the pressure gradient generated by the anti-cyclonic (less pro-grade) rotation. While the initial 
nebula (from 1 — 10 AU) had a mass of about IO^^Mq, there are only 8 x IO^^Mq left after 10^ 
years which corresponds to a mass loss (radially inward and outward) of 2.0 x 10“'^ Mq yr“^. We 
cannot tell in this simulation how fast mass is being accreted onto the star as our computational 
domain ends at 1 AU. The red blob collects a total of 10^° g (170 M 0 or 0.5 Jupiter masses 
[Mj]). Without further addition of mass or cooling this object is not gravitationally unstable, as 
the Toomre Q in the disk is about 10 and the local Jeans mass in the condensation is about 2.5 
Mj. 

We see that even in a disk which is not massive enough to fragment into planets or brown 
dwarfs (Boss 1998), a kind of pre-protoplanet can form simply as a result of baroclinicity and the 
resulting vorticity. The object, which is not yet gravitationally bound, could evolve into a planet in 
one of two ways: (1) it could efficiently collect solid particles in the center and wait until the critical 
mass for gas accretion is reached, or (2) it could concentrate enough gas and cool down efficiently 
to become gravitationally unstable. In either scenario the time scale for planet formation will be 
shorter than in cases without the vortices in the disk. Additionally the vortex scenario can explain 
why there could be a solid core in a planet even it formed by gravitational instability rather than 
by accretion of solid material to form a core followed by gas capture. 
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7. Conclusions 

The global baroclinic instability is found to generate turbulence in disks and drive an accretion 
process. 

In general we found numerically that isotropic turbulence in the R — (p plane of the disk has the 
property of transporting angular momentum radially outward. But only the global baroclinic insta¬ 
bility seems to provide a reliable source for this turbulence in the first place. Thermal convection 
in the vertical direction of the disk is not necessary for this effect. 

These results seem to be unaffected by the type of simulation. Whether 2D or 3D, whether 
polytropic, artificially heated or not heated, whether open boundaries or shearing disk, all simula¬ 
tions lead to the same conclusions. Nevertheless, only global 3D non-heated models will have the 
credibility to predict exact a-values. 

We showed that a protoplanetary disk with a density and temperature distribution which 
cannot be described by a single polytropic K is not barotropic and is possibly unstable against 
non-axisymmetric perturbations. This non-isentropic situation (with a radial entropy gradient) is 
ultimately a consequence of the radially decreasing gravitational force. The vertical component of 
gravity pushes the disk together and determines the pressure. Thus it is natural for non-isothermal 
disks to be baroclinic. 

We demonstrated numerically in a simple 2D simulation that a baroclinic disk is unstable and 
develops strong geostrophic turbulence while a barotropic disk is perfectly stable. 

Now we can answer the questions raised in §4. (A): Balbus et al. (1996; see also Hawley 
et al. 1999) could not observe this kind of instability as their simulations were barotropic. The 
simulations in SB96 allowed for local baroclinic effects but no global entropy gradient was present. 
They used the shearing-box approximation; thus f3p, (3 t = 0 and (5k is automatically zero. (B): 
The instability is purely hydrodynamical with an initial e-folding time for the growth of about 5 
orbital periods for an entropy gradient with (3k = 0.57. It occurs naturally in rotational shear 
flows when surfaces of constant density are inclined versus surfaces of constant pressure. A detailed 
stability analysis does not exist yet, but this is true for a lot of turbulent situations. Indications 
are, however, that the lowest-order modes have the fastest growth rates. It is also not known 
yet whether we are observing a linear or a non-linear instability operating in the disks. (C); The 
barotropic shearing-disk simulations as well as the baroclinic open-boundary simulations indicate 
that the “shearing disk” boundary conditions are not responsible for the turbulence. 

Thermal convection is only indirectly related to the baroclinic instability, as convective and 
radiative transport are responsible for the radial temperature distribution and therefore the radial 
entropy gradient. An isothermal disk would always be barotropic, but we are considering only the 
situation where the optical depth is larger than one, so the disks are not isothermal. The optical 
depth greater than one is also necessary for thermal convection, but a vertically convectively stable 
purely radiative disk can also establish a radial entropy gradient. If convection is present, then it 
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can be important in generating the initial non-axisymmetric perturbations induced in the Keplerian 
disk which are necessary to set off the baroclinic instability. On the other hand convection produces 
negative Reynolds stresses; however these turn out to be orders of magnitude weaker than the ones 
created by the baroclinic instability. 

The simulations also generate vorticity, as can be made plausible by a simple argument. Equa¬ 
tion (7) shows how a density fluctuation in the direction leads to a pressure gradient that does not 
exactly line up with the density gradient. Thus imagine two parcels of gas on the same orbit around 
the central object. When they get pushed apart in the azimuthal direction, then the pressure will 
try to restore the previous state by pushing them together again. In a barotropic disk they perform 
a damped oscillation, and the perturbation decays. But in a baroclinic disk, assuming that the 
perturbation occurred only along the azimuthal direction, the gas pressure will now not only push 
the parcels azimuthally together again, but also will drive the gas parcels slightly radially outward. 
For two reasons they are then driven back to the lower radius. First, there is the gas pressure at 
the larger radius, and, second, they do not have enough angular momentum to stay on that higher 
orbit. When they return to their lower orbit, they will fall behind their original azimuthal position, 
as the gas on a lower orbit has a larger angular velocity then the one on the outer orbit. And 
thus it is clear that a non-axisymmetric radial velocity distribution is created — dvrldcf) 7 ^ 0 — 
and thus vorticity generated. This leads first to a “meandering” flow which eventually becomes 
completely chaotic, characterized by irregular waves. Furthermore a radially local perturbation 
spreads quickly in the radial direction, as it always affects and perturbs neighboring radii. Other 
works (e.g. Adams &: Watkins 1995) have considered the behavior of vortices and their interaction 
with viscosity before, but none of them have created them. We show here for the first time that 
they form necessarily in a realistic disk in the absence of an underlying viscosity. 

Finally, the baroclinic instability generates turbulence with velocities close to the sound speed. 
Dissipation does not occur (homogeneously in time and space) on the small scales as in a models 
but occurs in large-scale shock structures. The calculations strongly suggest that this baroclinic 
instability is a feasible way to maintain turbulence and outward angular momentum transport in 
protoplanetary disks, even if the physical conditions do not allow for MHD turbulence. It would 
follow that there is no such thing as a “dead zone” in protoplanetary disks. Evidently then also the 
paradigm of layered accretion has to be revised. One also has to take into account a transition zone 
in a disk, which separates regions where MHD turbulence or hydro-turbulence dominate, and in 
which the two processes may coexist. It also would be fruitful and interesting to study the transition 
zone and in general the interaction between hydro turbulence and self gravitational instabilities, as 
such a process could be crucial for planet formation. These combined effects would be important 
at the earliest stages of evolution of a protoplanetary disk, when the disk is still relatively massive 
and material is still accreting onto the disk from the surrounding molecular cloud. Once the disk 
becomes optically thick and radiative transfer effects become important, one can expect baroclinic 
effects to occur. 

Gravitational and baroclinic instabilities seem to have certain properties in common. Each of 
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them leads to non-axisymmetric modes and an angular momentum transport that is not necessarily 
describable by an a— formalism. In baroclinic disks we are able to measure quite reasonable, and 
even more significant, all positive values for the Reynolds stresses, but it could be dangerous to use 
these values in the viscous description of an accretion disk. Our global model shows dramatically 
how far a real disk can depart from the idealized axisymmetric laminar disk that evolves in quasi¬ 
steady state on viscous time-scales. Thus, results from a—disk models can only reflect the longterm 
average properties of disks. 

The trans-sonic turbulence would be expected to be very influential on the passive dust con¬ 
taminant in generating collisions and mixing as well as concentrating the dust grains. Nevertheless 
the Mach number associated with mixing of angular momentum is almost two orders of magnitude 
less than that of the turbulence, which reflects the non-isotropic 2D character of the turbulence. 
The spectral density (see Fig. 4 and Fig. 18) indicates that there is also turbulence on the small 
scales, but only resolution studies can show how much energy really decays towards the Kolmogorov 
scale. In short, a baroclinic disk is much more turbulent than ‘viscous’. 

Our simple global 2D model already shows azimuthal density fluctuations by a factor of four. 
But most up-to-date models which try to interpret observational data use axisymmetric disks. Our 
global simulation could be a first step in showing how the spectral (radiation) energy distribution 
(SED) of a baroclinically unstable disk would look like in contrast to an axisymmetric one. 

The anti-cyclonic rotating gas parcels are vortices, that could be the precursors of planetary 
formation. They can be thought of as pre-protoplanets. In this connection, more realistic 3D 
models with radiation transport could allow for higher mass concentration, as the cooling can 
locally decrease the pressure. The planets could form either by concentration of dust in the centers 
of the vortices, as was suggested by Tanga et al. (1996) and Godon & Livio (1999b), or by sufficient 
gas accretion onto a vortex so that it undergoes gravitational collapse (Adams &: Watkins 1995). 
This second process would happen in a similar fashion to the model by Boss (1998) with the big 
difference that the vorticity takes care of the local mass enhancement even after the disk or even 
parts of the disk have become gravitationally stable. It furthermore is an easy way to explain a 
solid core for Jupiter and Saturn, as the pre-planetary embryo, the rotating slightly over-dense 
vortex, will already have accumulated planetesimals. 

In a final speculation we suggest that there is a connection between UXor events (Natta et al. 
1999) and the baroclinic instabilities. It would be logical to expect over-dense vortices to form in 
disks around Ae/Be-stars. In a nearly edge-on disk, a transit of such a cloud, which has a higher 
scale height than the material in the surrounding disk, could obscure the stellar light. In a viscous 
a—disk these over-dense regions would smear out in the azimuthal direction on a dynamical time- 
scale, which is an orbital period, but in a viscosity free disk, they do not only persist, they can 
form. 

Further 2D and 3D simulations are necessary to determine the proper role for the global 
baroclinic instability in accretion disk theory. Questions to be addressed are abundant, for example 
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what is the influence of the aspect ratio H/R oi the disk on the instability? What are the critical 
values for /3x? Can one measure the growth rates for idealized perturbations with a single wave 
number m? 

Apparently the dissipation and conservation properties of different numerical hydro schemes 
affect the development of the instability. We already noticed that not all available codes show the 
same results if one does baroclinic simulations. We will continue with these tests and suggest that 
other owners of hydro codes should try to do so as well, in order to prove or disprove our findings. 

The next step in our investigations of this new instability will be a detailed stability analysis, 
which has to be the ultimate proof for our numerical findings. 
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name 

grid 


S(p 

Pk 

P 

II 

M 

^intf tdyn 

Model 1 

51 X 20 X 60 

3.5~6.5 

90° 

- 

2 X 10-3 

0.5 

68 

Model IB 

51 X 20 X 60 

3.5-6.5 

90° 

- 

2 X 10-2 

1.5 

14 

Model 2 

64^ 

o 

1 

o 

30° 

0.0 

2 X 10-^ 

10-^ 

no 

Model 2B 

64^ 

o 

1 

o 

30° 

0.0 

3 X 10-^ 

0.04 

1 

Model 3 

64^ 

o 

i 

o 

30° 

0.57 

1.5 X 10-^ 

0.03 

no 

Model 4 

1282 

o 

1 

o 

30° 

0.57 

2 X 10-^ 

0.04 

130 

Model 5 

1282 

CO 

o 

i 

o 

60° 

0.57 

2.5 X 10-^ 

0.05 

63 

Model 6 

1282 

1.0-10.0 

360° 

0.57 

8 X 10-3 

0.3 

50 


Table 1: Parameters chosen in the different simulations. These parameters are the dimensioning of 
the grid {nr,n^,nip) or {nr,n^), the radial extent of the disk Ri to Rq-, the azimuthal extent 5ip, 
the parameter (Sk which describes the variation of the polytropic K with radius, the strength of 
the measured Reynolds stresses a = Ar^j), the Mach number (M) of the turbulence, and the ratio 
of the time over which the averaging was performed to the orbital period at the outer edge of the 
disk. 



- 28 - 


Fig. 1.— Three-dimensional rendering of the vertical mass flux in model 1, at a time near the end 
of the run. Red denotes positive velocities, and blue denotes negative. Compare this figure to Fig. 
3 in Stone & Balbus (1996). 

Fig. 2.— Turbulence in model 1: Vertically-, azimuthally- and time-averaged stresses over 68 orbits, 
measured in units of an effective A (see eq.[ 22]) and plotted as a function of radius (in AU). The 
mean value for the angular momentum transport {upper left: a = Arcf, ~ 2 x 10“^) is given by the 
dotted line. Other frames give the averages of the relative density fluctuation, the strength of the 
turbulence in terms of velocity fluctuations (eqs.[24], [25], and [26]), and Mach number. See §3.2 
for explanation. 

Fig. 3.— Model 1 at a time near the end of the run; Surface density ( colors; 96 [violet] to 335 
[red] g/cm^), velocities (vectors; Vmax = 0.6x the local sound speed) and iso-temperature contours 
in the midplane. 

Fig. 4.— Model 1 at the same time as in Fig. 3; Spectral density distribution of the velocities at 
the midplane computed along the (/^-direction and averaged over radius. The slope for isotropic, 
incompressible turbulence (i.e. a Kolmogorov spectrum) is indicated by the dashed line and the 
spectrum for 2D geostrophic flows by the dotted line. 

Fig. 5.— Model IB; Surface density 52 orbits after the initial state ( colors; 20 [violet] to 523 [red] 
g/cm^), velocities (vectors; Vmax = 0.75x the local sound speed) and iso-temperature contours in 
the midplane. 

Fig. 6.— Turbulence in model IB; Vertically-, azimuthally- and time-averaged stresses over 14 
orbits, taken 52 orbits after the initial perturbation, measured in units of an effective A (see eq.[ 
22]) and plotted as a function of radius (in AU). The mean value for the angular momentum 
transport {upper left: a = w 2 x 10“^) is given by the dotted line. Other frames give the 

averages of the relative density fluctuation, the strength of the turbulence in terms of velocity 
fluctuations (eqs.[24], [25], and [26]), and Mach number. See §3.2 for explanation. 

Fig. 7.— Model 1; Radial distribution of the polytropic K (normalized to the biggest value) as 
measured from the temperature and density in the midplane in the thermal convective flow before 
the system became turbulent, i.e. it was axisymmetric {dotted lines: snapshots at 10 different times; 
solid line: theoretical value for (3k = 0.57). 

Fig. 8.— Model 2; Surface density (colors; 277 [violet], 300 [green] to 305 [red] g/cm^ [same color 
coding as in Fig. 12]) after 90 orbits, velocities (vectors; Vmax = 1 x 10“^x sound speed) and iso¬ 
temperature contours in the midplane. Barotropic simulation with no growing non-axisymmetric 
instability. 

Fig. 9.— Model 2; Time development of kinetic energy in the r-direction {upper left), in the (j)- 
direction {upper right), spatially mean accretion rate in Mq yr“^ averaged over five orbits {lower 
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left)^ and enstrophy (integral square vorticity; lower right). The units for kinetic energy and 
enstrophy are normalized to the first value occurring in the data set. As predicted for barotropic 
flows, no instability growth can be observed within the first 100 orbits. No vorticity is generated. 

Fig. 10.— Model 2: Azimuthally- and time-averaged Reynolds stress over 10 orbits {upper left) 
taken 100 orbits after the initial perturbation, measured in units of an effective Ar^ (see eq. [22]) 
and plotted as a function of radius (in AU). The mean value is a = Ar^f, ~ 2.0 x 10“®, which is 
practically zero. Other frames give the overall Mach number, and the strength of the turbulence 
in the radial direction (Arr; see eq. [24]) and the azimuthal direction {A^^; see eq. [26]). 

Fig. 11.— Model 2B: Azimuthally- and time-averaged Reynolds stress over 1 orbit {upper left) 
right after a strong initial density perturbation, measured in units of an effective (see eq. [22]) 
and plotted as a function of radius (in AU). The mean value is a = A^^ ~ 3 x 10“^. Other frames 
give the overall Mach number, and the strength of the turbulence in the radial direction {Arr', see 
eq. [24]) and the azimuthal direction (A,^^; see eq. [26]). 

Fig. 12.— Model 3: Surface density (colors; 286 [violet] to 318 [red] g/cm^), velocities (vectors: 
Vmax = 0.03 X sound speed) and iso-temperature contours in the midplane after 600 orbits. Baro- 
clinic simulation with a still growing non-axisymmetric instability. 

Fig. 13.— Model 3: The quantities plotted have the same meaning as those in Fig. 9 and are 
directly comparable. In baroclinic flows, vorticity and enstrophy are created and instability growth 
can be observed within the first 100 orbits. 

Fig. 14.— Model 3; Azimuthally- and time-averaged Reynolds stress over 380 orbits (starting with 
2030 orbits after the initial perturbation), measured in units of an effective Artp (see eq. [22]) and 
plotted as a function of radius (in AU). The mean value {upper left) is a = A^^ 1.5 x 10“*^. Other 

frames show the overall Mach number and the strength of the turbulence in the radial direction 
{Arr', see eq. [24]) and in the azimuthal direction (A^^; see eq. [26]). 

Fig. 15.— Model 4: Stresses averaged over 380 orbits (starting with 2030 orbits after the initial 
perturbation). The quantities plotted have the same meaning as in Fig. 14. The mean value of 
a = Ar^ ^2 X 10“^. 

Fig. 16.— Model 5: Stresses averaged over 180 orbits (starting with 600 orbits after the initial 
perturbation). The quantities plotted have the same meaning as in Fig. 14. The mean value of 
a = Ar^ ~ 5.0 X 10“^. 

Fig. 17.— Model 5: The quantities plotted have the same meaning as those in Fig. 12. This 
calculation was run over 230 orbits. 

Fig. 18.— Model 5: Spectral density distribution of the velocities at the midplane computed along 
the (/9-direction and averaged over radius. The slope for isotropic, incompressible turbulence (i.e. a 
Kolmogorov spectrum) is indicated by the dashed line and the spectrum for 2D geostrophic flows 
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by the dotted line. 

Fig. 19.— Model 6: Four snapshots of the evolution of surface density (colors: 650 [red], 550 
[yellow], 450 [green], 250 [blue] to < 100 [black] to g/cm^) in the global model in polar (r — (p) 
coordinates. The times are a): 1, b): 95, c); 190 and d); 320 orbits at the outer radius. Contours 
of equal pressure are also shown. 

Fig. 20.— Model 6; Stresses averaged over 50 orbits (starting with 190 orbits after the initial 
perturbation). Quantities plotted have the same meaning as in Fig. 10. The mean value of Aj.(p ~ 
8.0 X 10"3. 

Fig. 21.— The “pre-protoplanet” in Model 6: Surface density (colors: 650 [red], 550 [yellow], 450 
[green], 250 [blue] to < 100 [black] to g/cm^) in the global model is projected in a cartesian frame 
after 320 orbits at the outer radius, which corresponds to 10^ yr. Note that the condensation is 
partially artificially smeared out in the (/^-direction, which is a result of the low order advection 
scheme. In reality one could expect the “pre-protoplanet” to be more strongly confined. 

Fig. 22.— Model 6: Surface density (colors: 650 [red], 550 [yellow], 450 [green], 250 [blue] to 
< 100 [black] to g/cm^) and velocity (vectors: Vmax ~ l-5x sound speed) in the global model in 
polar (r — tp) coordinates after 320 orbits at the outer radius, which corresponds to 10^ yr. Plotted 
velocities in the p direction are obtained by subtracting the mean azimuthal velocity at each radius, 
which explains why the vortex center seems to be displaced from the density maximum. 
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